Linear response formula for open systems 
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An exact expression for the finite frequency response of open classical systems coupled to reservoirs 
is obtained. The result is valid for any conserved current. No assumption is made about the 
reservoirs apart from thermodynamic equilibrium. At non-zero frequencies, the expression involves 
correlation functions of boundary currents and cannot be put in the standard Green-Kubo form 
involving currents inside the system. 

PACS numbers: 



I. INTRODUCTION 

One of the important tools in the study of transport 
phenomena is the Green-Kubo formula[H 0|, which re- 
lates the equilibrium correlation function of any two con- 
served currents (i.e. currents associated with conserved 
charges) to the linear response of the first when a gradi- 
ent in the potential conjugate to the second is applied. 
The response is expressed in terms of conductivity co- 
efficients of a system, in the thermodynamic limit when 
boundary effects can be neg lected. The proof@ relies on 
Onsager's[3| relation between the time evolution of equi- 
librium and non-equilibrium fluctuations, or obtains the 
response to an external field [3 and relates it to the re- 
sponse to an internal gradient [5[ • One can thence obtain 
the conductance of a system if it is sufficiently large and 
if the conductivity does not diverge. 

If either of these conditions breaks down, one has to 
find the conductance of the system directly. This can 
be done by connecting it to infinite reservoirs [f| 0], but 
requires assumptions about the reservoirs and the appro- 
priate ones can be subtle @. Alternatively, a Green-Kubo 
like formula has been obtainedfl] for the conductance(s) 
of an open classical system without having to deal with 
the reservoirs. This is done by assuming that the system 
is in a non-equilibrium steady state and (implicitly) that 
the current driven by the external potential gradient is a 
scalar. The assumptions restrict the result to the zero- 
frequency conductance of a quasi one-dimensional (in a 
sense to be made precise later) system. 

In this paper, we derive a formula for the generalized 
linear conductance of a finite classical system without any 
of these restrictions. The system can have an arbitrary 
shape and number of reservoirs, and the formula applies 
at any frequency. The only assumption made about the 
reservoirs is that they are in thermodynamic equilibrium. 
The coupling between the system and the reservoirs is 
assumed to be such that if the system starts in equilib- 
rium at the same thermodynamic potentials as all the 
reservoirs, it remains in equilibrium, but we do not as- 
sume that the reservoirs can equilibrate the system. For a 
quasi one-dimensional system at zero frequency, the for- 
mula can be transformed into that of Ref. [9( which is of 
the standard Green-Kubo form, i.e. it involves the equi- 
librium fluctuations of currents inside the system. We 



show that this transformation is not possible at finite 
frequencies. 

The derivation here builds on previous work that ob- 
tained the thermal conductance at zero .1 Of and non- 
zero [llj frequency of a finite classical system in con- 
tact with heat baths. Although a large variety of heat 
baths were considered, the proof had to be painstakingly 
constructed separately for each bath. Thus it was not 
clear whether it might fail for some types of heat baths 
and whether it applied to other conserved currents. The 
present paper resolves these questions. 

The rest of this paper is organized as follows. In Sec- 
tion fll] , we briefly review the result of Ref. [T3, [H| f° r 
the thermal conductance of a Fermi-Pasta-Ulam (FPU) 
chain with Langevin heat baths at the ends. The no- 
tation is more general and part of the proof is slightly 
different to allow it to be extended to other conserved 
currents. In Section [TTTl we extend the proof to the case 
of all conserved currents. In Section HVl we transform the 
formula at zero frequencies to the 'standard' form, and 
discuss why this is not possible at non-zero frequencies. 



II. FPU CHAINS WITH LANGEVIN BATHS 

We first review the derivation of Refs. [l(| El for 
the heat conductance of a iV-particle FPU chain with 
Langevin baths at the ends. The equations of motion are 

d 

mik = —3— [U{xi-x - x{) + U(xi - xi+i) +V(xi)\ 
0x1 

+ Si,i[r)L(t) - jlvi] + 6i,N[r)R(t) — ~/rvn] (1) 

for I = 1, 2, ... N. Here to/, xi, vi are the mass, position 
and velocity of the Z'th particle. U and V are the inter- 
particle and onsite potentials with xq = xn+i = 0. jl,r 
and r/L,R(t) are the damping and Gaussian noise from 
reservoirs at temperatures Tl.r, satisfying 



(2) 



In the first stage of the proof, the Fokker Planck 
equation for the full phase space distribution function 
P(x;v;i) is constructed, where x = {xi...xn} and 
v = {vi . . . vn}. If Tl = Tr,, the steady state solution to 
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the equation is the equilibrium Boltzmann distribution 
P°(x, v). For T L<R = T± ±AT, we have 



dP 
~dt 



LP + L AT P 



(3) 



where L is the equilibrium (AT = 0) Fokker Planck op- 
erator 
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mi dv\ dvi 
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+ (l,i) -> (N,R) 



d 

1 dx; 
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(Fi is the force on the i'th particle) and 



-AT _ k B AT 



1L— _ lR d2 

m 2 dvf m 2 N dvj 



(4) 



(5) 



With T(x, v, t) = P° + p(x, v, t), to linear order in AT 



p(x;v;t)=/ e (t -^ L r fp (v)P°( X ,v)Ap(t')dt' (6) 



where A/3 = A(l/fc B T) = -AT/{k B T 2 ) and J e fp {w) is 
defined by 



<9P 
~dt 



L AT P° = (A/3)J? p P° 



(J) 



p=p° 



(the superscript in Jj referring to the energy). For any 
observable A, we define (5 A) = (A) — (A)o in terms of its 
expectation values with Af3 ^ and A/3 = 0. Then 



(SA(ti)) AT = / (A(ti)J%{t 2 )) eq Ap{t 2 )dt 2 (8) 

where the correlation function on the right hand side is 
evaluated in equilibrium with AT = 0. 

More generally, if the reservoirs at the ends have differ- 
ent values for some thermodynamic potential $ p whose 
conjugate conserved density is p, then 



(5A(ti)) 



A* 



(A(ti)J p fp (t 2 )) eq A<pp(t 2 )dt 2 (9) 



with Jj p defined by the generalization of Eqs.Q and ([7]). 
For example, for the particle current, $™ = — /3/i, where 
p, is the chemical potential. 

The second part of the proof is specific to Langevin 
baths at different temperatures, and is slightly different 
from Refs. [10|, 111) . From the Fokker Planck equation, 
one can verify that 



fp 



lR 
2mN 



[m N v% - k B T] 



1l 
2mi 



[miv 2 -k B T]. (10) 



We define the boundary energy current variable as the 
mean of the instantaneous energy currents flowing into 



the system from the left reservoir and flowing out of the 
system to the right reservoir. Thus 



(11) 



where 



(12) 



fixit) = -Ti«i 2 (i)+^(tW(«), 

i^,Ji(*) = -lRV 2 N {t)+TlR{t)v N {t). 

Then we will prove the relation 

(A(ti)J e fp (t 2 )) eq = -(A(ti)Ji(t 2 )) eq if h > t 2 . 

(13) 

To prove Eq. fp~3|) . we use Eq. ([l~2"T) and average over the 
noise. In performing this average, we discretize the equa- 
tions of motion Eq.([T]) in the usual manner, with the right 
hand side evaluated at time t and the left hand side equal 
to mi[vi(t + e) — vi(t)]/e. The boundary heat currents are 



— 1 

— 2 



Ui(t)+«i(t + e)][»fr(t)-7iUi(t)] 
[v N (t) + v N (t + e)][r) R (t) ~ !RVN(t)]{U) 



Usually, the noise average is performed at fixed v\{t) 1 
yielding 



(3lL(t))v=lL(k B T/mi-vl(t)). 



(15) 



Naively, this would seem to yield Eq. ([T3| without the 
minus sign on the right hand side. 

However, we want to use the Fokker-Planck evolution 
operator over the interval ti > t > t 2 to evaluate the 
left hand side of Eq. (fT3"l) . Therefore, we use J^(t 2 — e) in 
Eq. ([13[) and perform the noise average for fixed V\ jv(*2), 
not vi } N{t 2 — e). We change variables in Eq. (fT4]) from 
VL,R(t) to fj LtR (t) = rj LjR (t) - 2jl,rVi,n, since r] LM (t) 
are uncorrelated with Ui,jv(t + e). It is easy to verify that 



*))v 



lL{v{{t)- 
lR^ 2 N {t) 



k B T/mi), 
- k B T/m N ) 



(16) 



The same result can also be obtained by using the fact 
that evolving backwards in time from t 2 reverses the cur- 
rents, and from the time reversal invariance of the dy- 
namics the noise-averaged time reversed current is given 
by Eq. (fT5)) . Comparing with Eq.©, the noise averaged 
Jb(t2 — e) is equal to —Jf(t 2 ), thus proving Eq. (fT3)) . 
Combining Eqs.© and (fHf]) . we have 



(SA(ti)) 



AT 



<,4(ti)J e &(t 2 )) B ,A/3(t 2 )dt 2 (17) 



which gives the response to a general time dependent (i.e. 
non-zero frequency) variation in the temperatures of the 
reservoirs. Equivalently, Fourier transforming, 



(SA(u))at = -A/3(w) 



ex.p[iut](A(t)J§(0)) eq dt 



This equation is now generalized to all conserved currents 
in Section Hm 
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III. GENERALIZED CURRENTS 



side of Eq. ([20|) can then be transferred: 



We first consider the case of particle currents. Al- 
though Section [TT| was for a lattice system with a fixed 
number of particles, it can be easily extended to a sys- 
tem with a continuous coordinate x in which the particles 
move around, allowing particles to be enter and leave the 
system from the reservoirs. 

Because the number of particles in the system is no 
longer conserved, we work in the grand canonical ensem- 
ble. Let P/v(x, v, t) be the probability density for the sys- 
tem to be in an iV-particle configuration with coordinates 
(x, v) at time t. Thus J Piv(x, v, t)dxdv is the probability 
for the system to be in an TV-particle configuration. In 
equilibrium, J P$(x, v)dxdv = exp[/3/xiV]Zjv/0, where 
Zn is the canonical partition function for the A^-particle 
system and Q is the grand partition function. The Fokker 
Planck probability density is now an infinite column vec- 
tor, and the time evolution operators L and L Ap are 
matrices. We will denote the column vector as P(x, v, t). 
Analogous to Eq.©, we define JJ through 



dP 
~dt 



P=P» 



(19) 



The operator L AfJ - is due to the difference in the chemical 
potentials of the reservoirs and the system, and changes 
the number of particles N. 

Let n (/3[il) be the transition rate from 

Pn(xi ■ ■ ■ v n) to Ppfi (x[ . . . v' N i) due to the left reservoir. 
(The dependence of A on other variables such as x, v 
and t is implicit; if some of the reservoir variables are 
so slow that they cannot be integrated out to yield 
transition rates, one can augment the arguments of P/v 
and the A's to include these variables.) The contribution 
from the left reservoir to the first order equation for 
APv = P N -P% is 

^AP ff ( Xl ...v N ) = A(/3 M ) L ]T f [(d^ L X L NrN ,)P^ 

N' J 

- (^ li A^)P°]dx' 1 ...dvV(20) 

where the derivatives on the right hand side are taken 
at (/3/i)l = Pn- Here <9 t L AP/v is the part of the time 
evolution of AP/v due to the left reservoir; d t PN = (dt + 
d*)P N . 

Since the left reservoir cannot disturb the grand canon- 
ical ensemble distribution when it is at the same j3fj, as 
the system, we have 



E 

N' 



\ L P° /fir' 

A N,N'- r N' ax l 



.dvv = E/ 

N' J 



W' N 



P^dxi ...dv'j 



N 1 



whenever (/3/i)z, = fil*- The derivatives on the right hand 



5 t i APv(x 1 ...v Ar ) = -A(/3fi) L ^2 I 

N' 



dx[ . . . dv' N (22) 



Since P^ ~ exp[iV/3/i], applying Eq. ([2T|) again 



d^AP N = A(j3^) L J2(N - N') 



N' 



\L p 
A N.N'- r N' 



(23) 



Now A^r n'Pn' integrated over x^ . . . v' N , is the rate at 
which systems transition from an N' particle state to 
an A-particle state at the phase space point xi...vjv- 
Therefore the probability density for the system to be 
in an A-particle state at time t and have received N — 
N' particles from the left reservoir between time t — St 
and t is N ,P^,8t. For a system in an A-particle state 
at Xi . . . v^r, the particle current from the left reservoir 
immediately before the time t is then j^ v i P ( ^ v (x, v) = 
J2n'(N — N')\j^ N /Ppf'- The stipulation that the current 
has to be evaluated just before t is important as seen 
from the discussion around Eqs. (TT5|) and (TTB)) . 

Comparing with Eq. ([23|) . including the effect of the 
reservoir to the right with A(/3/x)l = —A(f3^)n = 
\A(Pfj) and using Eq. ([T9j) . we obtain 



(24) 



where J£ is defined in terms of the particle currents from 
the reservoirs in the same manner as Eq. (fTTj) . 

Although energy currents were handled differently in 
Section [U there is no reason why one could not have 
proceeded as we have done here: define Pe(x, v, t) as the 
probability density for the system to have an energy E 
and be at the phase-space point (x, v) . (The sums over N' 
would have been replaced by integrals over E 1 .) Since the 
number of variables in (x, v) is independent of the energy 
E and different energies correspond to non-overlapping 
regions in phase space, this would have been an unnec- 
essary complication. But from this we see that although 
the specific example of particle currents has been used in 
this Section, all conserved currents can be treated in the 
same manner, and 

(A{h)J p fp {t 2 )) eq = -(A(ti)Jg(t2)U if h > t 2 (25) 

for any conserved current J p . Nor is one limited to 
two reservoirs: since the proof we have constructed for 
Eq. (|25j) deals with the reservoirs separately, one can con- 
sider an arbitrary number of reservoirs at potentials <£>£ . 
Combining with Eq.©, 

(6A(t)) A *=J2 f (A{t)f ba {t')) eq AK{t')dt' (26) 

„, J — OO 
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If A is the current j% at the boundary k, one obtains 
the conductance 1131: 



O b %W>A* = E f {3lAt)f b , a { tr ))^K{t')dt'. 

a,p J -°° 



(27) 



As with Eq. ffTT)) this can be Fourier transformed to yield 

poo 

0'6,»>A* = J2 A$ °M / °fpM{jM(04(o))^ 

(28) 

The equilibrium correlation function in Eq. (|26p must be 
evaluated with the same boundary conditions as the left 
hand side, i.e. with reservoirs connected to the open sys- 
tem (but, unlike the left hand side, with the reservoirs at 
the same thermodynamic potentials). Different bound- 
ary conditions may be inequivalent even in the large sys- 
tem limit [Til [l3lj. 



IV. ZERO FREQUENCY CONDUCTANCE 

At zero frequencies, for a quasi one dimensional system 
with two reservoirs, Eq.Q can be transformed into a 
more familiar form. The proof is a generalization of the 
one in Ref. [l(j. If the reservoirs are at potentials $ p ± 
(A$ p )/2, from Eqs.© and (|25]) 



(SA(t)) 



A* 



P 



(A(t)J p (t')) eq dt'. (29) 



For concreteness, we first consider a one dimensional 
lattice. Let pi be the amount of the conserved quantity 
p at the i'th particle, i be its current between the 

i'th and i + l'th particles, and J p {t) — YliLo^i+i ^ 
A = J CT is the current flowing through the entire chain 
that is associated with a conserved quantity er, Eq. (|2"9"|) 
becomes 



(8J°{t)) 



A* 



E 

p 



(A$ p 



(r(t)jp(t')) eq dt'. (30) 



Define D p (t) = J2k=i( N + 1 - 2fc)p fe . From the conti- 
nuity equation, 



dD p /dt = -2J p (t) + 2{N - 1) J h p (i). 



(31) 



Multiply both sides of this equation by J a (ti), take an 
equilibrium average, and integrate over — oo < t < t±. 
Then {J°{t)DP {-oo)) eq = (J a ) eq (D p ) eq = 0, and if p 
and a have the same symmetry under time reversal or the 
system is symmetric under reflection about its middle, 

{J a (t)D p {t)) eq =0. (32) 

If either condition is satisfied for each p in Eq. ([3U)) , 

1 



(5J°{t)) 



A* 



N — 1 ^ 

p 



(A$") / {J°{t)JP(t')) eq dt'. 



a=l,2 



■P 



For a continuum system of length L, D p is defined 
as J p(x)(L — 2x)dx, and Eqs. ([3T|) and (133)) are obtained 
with N-l replaced by L. The factor of A — 1 in Eqs. ([3"Tj) 
and p3p is replaced with L. In higher dimensions, we 
define D p = —2 J xp(r)dr and use 

dD p /dt = 2 J xV- j p (r)dr = -2 J j p {r) + 2 

(34) 

where ja; is the x component of the current and the sum is 
over the reservoirs. Then if Eg. (f3"2"|) is satisfied, Eg. (1531 is 
obtained with J p — > J p on the right hand side and N — 1 
replaced by the x-separation of the two reservoirs [lij. 

Eg. (|32p is always satisfied for quasi one dimensional 
homogeneous systems. By quasi one dimensional, we 
mean that the system is a tube whose cross-section does 
not vary along its length, with reservoirs at its ends. 
Translational invariance within the tube ensures that cur- 
rents in response to potential differences between the 
reservoirs are parallel to the orientation of the tube, 
and so are effectively scalar. This is implicitly assumed 
in Ref. [9(, which derives Eq. (|33[) for a non-equilibrium 
steady state. Eq. ([32l) is also satisfied if p = a. 

Since both particle and energy currents change sign 
under time reversal, the continuum version of Eq. (|33[) is 
valid for thermoelectric transport coefficients: 



f = G nn A(/3p) - G ne A(3 
f = G en A(l3p) - G ee A(3 



(35) 



(33) 



where j , j are non-equilibrium steady state currents 
flowing in response to differences in the chemical po- 
tentials and temperatures of the reservoirs, and G ap = 

f™{r(t)j p (o))dt. 

Despite the similarity between Eq. (13"3"|) and the Kubo 
formula for the conductivity @, H| , there are important 
differences. The caveat at the end of Section Mil is still 
valid: the correlation function has to be calculated with 
the correct boundary conditions. Second, Eg. ([33)1 has 
only been obtained in the zero frequency limit. 

One might wonder whether the w ^ version of 
Eg. ([33)) . with A$ p (i') taken inside the integral, might be 
true even if not proved here. A simple example proves 
otherwise: if the contacts to the reservoirs are so weak 
that the system is effectively isolated, a non-zero A$ 
produces no response and the boundary current j£ in 
Ea. ([21))) is zero. However, (J CT (t)J p (t') eq ) ^ because 
currents flow inside an isolated system due to sponta- 
neous fluctuations, and ^_ oo (J a {t)J p (t') eq A^ p {t')dt l ^ 
0. By contrast, if the system is driven by an applied field 
instead of a difference in potentials between the reser- 
voirs, the w / version of Eg. ([33)) is always valid. 

To summarize, we have obtained a formula for the gen- 
eralized finite-frequency response of a finite classical sys- 
tem connected to an arbitrary number of reservoirs when 
the thermodynamic potentials of the reservoirs are var- 
ied. The equilibrium correlation function in the formula 
involves currents flowing in from the reservoirs. At zero 
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frequency, it reduces to a familiar Green-Kubo form, i.e. the system but this is not the case at finite frequencies, 
it involves equilibrium fluctuations of the currents inside 
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